Step 1: Define the problem

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn')

import warnings
warnings.filterwarnings('ignore')
warnings.filterwarnings(action='ignore', category=DeprecationWarning)
warnings.filterwarnings(action='ignore', category=FutureWarning)

missing_values = ["?", ".", "", "_", "Na", "NULL", "null", "not", "Not", "NaN", "NA", "??", "nan", "inf"]
raw_data = pd.read_csv("Police_Department_Incident_Reports__2018_to_Present.csv", na_values=missing_values)
raw_data.set_index('Row ID', inplace=True)
raw_data.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 410210 entries, 95308704134 to 95316305151
Data columns (total 25 columns):
 #   Column                   Non-Null Count   Dtype  
---  ------                   --------------   -----  
 0   Incident Datetime        410210 non-null  object 
 1   Incident Date            410210 non-null  object 
 2   Incident Time            410210 non-null  object 
 3   Incident Year            410210 non-null  int64  
 4   Incident Day of Week     410210 non-null  object 
 5   Report Datetime          410210 non-null  object 
 6   Incident ID              410210 non-null  int64  
 7   Incident Number          410210 non-null  int64  
 8   CAD Number               318808 non-null  float64
 9   Report Type Code         410210 non-null  object 
 10  Report Type Description  410210 non-null  object 
 11  Filed Online             84136 non-null   object 
 12  Incident Code            410210 non-null  int64  
 13  Incident Category        409921 non-null  object 
 14  Incident Subcategory     409921 non-null  object 
 15  Incident Description     410210 non-null  object 
 16  Resolution               410210 non-null  object 
 17  Intersection             389130 non-null  object 
 18  CNN                      389130 non-null  float64
 19  Police District          410210 non-null  object 
 20  Analysis Neighborhood    389043 non-null  object 
 21  Supervisor District      389130 non-null  float64
 22  Latitude                 389130 non-null  float64
 23  Longitude                389130 non-null  float64
 24  point                    389130 non-null  object 
dtypes: float64(5), int64(4), object(16)
memory usage: 81.4+ MB

raw_data.head()

	Incident Datetime	Incident Date	Incident Time	Incident Year	Incident Day of Week	Report Datetime	Incident ID	Incident Number	CAD Number	Report Type Code	...	Incident Description	Resolution	Intersection	CNN	Police District	Analysis Neighborhood	Supervisor District	Latitude	Longitude	point
Row ID
95308704134	15/08/2020 12:43	15/08/2020	12:43	2020	Saturday	15/08/2020 12:58	953087	200490354	202281583.0	II	...	Battery	Open or Active	GENEVA AVE \ LONDON ST	21475000.0	Ingleside	Excelsior	11.0	37.716039	-122.440255	(37.716038818883085, -122.44025513581519)
64999771000	18/01/2018 19:00	18/01/2018	19:00	2018	Thursday	22/01/2018 16:59	649997	186068683	NaN	II	...	Lost Property	Open or Active	NaN	NaN	Out of SF	NaN	NaN	NaN	NaN	NaN
95319604083	16/08/2020 3:13	16/08/2020	3:13	2020	Sunday	16/08/2020 3:14	953196	200491669	202290313.0	II	...	Firearm, Discharging in Grossly Negligent Manner	Open or Active	23RD ST \ ARKANSAS ST	23642000.0	Bayview	Potrero Hill	10.0	37.754827	-122.397729	(37.75482657770952, -122.39772873392515)
95326228100	16/08/2020 3:38	16/08/2020	3:38	2020	Sunday	16/08/2020 4:56	953262	200491738	202290404.0	II	...	Malicious Mischief, Breaking Windows	Open or Active	VALENCIA ST \ 15TH ST	24377000.0	Mission	Mission	9.0	37.766540	-122.422044	(37.76653957529556, -122.42204381448558)
95322706244	15/08/2020 9:40	15/08/2020	9:40	2020	Saturday	15/08/2020 18:21	953227	206121692	NaN	II	...	Theft, From Locked Vehicle, >$950	Open or Active	NaN	NaN	Park	NaN	NaN	NaN	NaN	NaN

5 rows × 25 columns

raw_data.describe()

	Incident Year	Incident ID	Incident Number	CAD Number	Incident Code	CNN	Supervisor District	Latitude	Longitude
count	410210.000000	410210.000000	4.102100e+05	3.188080e+05	410210.000000	3.891300e+05	389130.000000	389130.000000	389130.000000
mean	2018.886119	802774.974791	1.904148e+08	1.912847e+08	25089.386536	2.532681e+07	5.965577	37.769380	-122.423775
std	0.789394	104865.461463	8.831851e+06	1.896881e+07	25785.005878	3.085594e+06	2.788198	0.024062	0.026120
min	2018.000000	618687.000000	0.000000e+00	1.000000e+00	1000.000000	2.001300e+07	1.000000	37.707988	-122.511295
25%	2018.000000	712073.250000	1.808363e+08	1.825131e+08	6244.000000	2.397300e+07	3.000000	37.756167	-122.434062
50%	2019.000000	802931.500000	1.904324e+08	1.914336e+08	7046.000000	2.491800e+07	6.000000	37.775873	-122.417707
75%	2020.000000	893632.000000	2.000466e+08	2.002821e+08	61030.000000	2.642200e+07	8.000000	37.785829	-122.407337
max	2020.000000	984190.000000	9.811720e+08	1.000000e+09	75030.000000	5.412200e+07	11.000000	37.829991	-122.363743

Step 2: Data cleaning / preprocessing & Feature Engineering

Loại bỏ các biến thừa không cung cấp thông tin hữu ích nào cho việc phân tích

data = raw_data.drop([
    'Incident Date'       , 'Incident Time'      , 'Incident Year'          , 
    'Report Datetime'     , 'Report Type Code'   , 'Report Type Description',
    'Incident ID'         , 'Incident Number'    , 'CAD Number'             , 
    'Filed Online'        , 'Incident Code'      , 'Incident Subcategory'   , 
    'Incident Description', 'CNN'                , 'Supervisor District'    , 
    'point'
], axis=1)
data.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 410210 entries, 95308704134 to 95316305151
Data columns (total 9 columns):
 #   Column                 Non-Null Count   Dtype  
---  ------                 --------------   -----  
 0   Incident Datetime      410210 non-null  object 
 1   Incident Day of Week   410210 non-null  object 
 2   Incident Category      409921 non-null  object 
 3   Resolution             410210 non-null  object 
 4   Intersection           389130 non-null  object 
 5   Police District        410210 non-null  object 
 6   Analysis Neighborhood  389043 non-null  object 
 7   Latitude               389130 non-null  float64
 8   Longitude              389130 non-null  float64
dtypes: float64(2), object(7)
memory usage: 31.3+ MB

Data Formating

data['Incident Datetime'] = pd.to_datetime(data['Incident Datetime'], format='%d/%m/%Y %H:%M')
data.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 410210 entries, 95308704134 to 95316305151
Data columns (total 9 columns):
 #   Column                 Non-Null Count   Dtype         
---  ------                 --------------   -----         
 0   Incident Datetime      410210 non-null  datetime64[ns]
 1   Incident Day of Week   410210 non-null  object        
 2   Incident Category      409921 non-null  object        
 3   Resolution             410210 non-null  object        
 4   Intersection           389130 non-null  object        
 5   Police District        410210 non-null  object        
 6   Analysis Neighborhood  389043 non-null  object        
 7   Latitude               389130 non-null  float64       
 8   Longitude              389130 non-null  float64       
dtypes: datetime64[ns](1), float64(2), object(6)
memory usage: 31.3+ MB

Loại bỏ các giá trị bị khuyết

data.isnull().sum().sort_values(ascending=False)

Analysis Neighborhood    21167
Longitude                21080
Latitude                 21080
Intersection             21080
Incident Category          289
Police District              0
Resolution                   0
Incident Day of Week         0
Incident Datetime            0
dtype: int64

data.dropna(axis=0, inplace=True)
data.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 388756 entries, 95308704134 to 95316305151
Data columns (total 9 columns):
 #   Column                 Non-Null Count   Dtype         
---  ------                 --------------   -----         
 0   Incident Datetime      388756 non-null  datetime64[ns]
 1   Incident Day of Week   388756 non-null  object        
 2   Incident Category      388756 non-null  object        
 3   Resolution             388756 non-null  object        
 4   Intersection           388756 non-null  object        
 5   Police District        388756 non-null  object        
 6   Analysis Neighborhood  388756 non-null  object        
 7   Latitude               388756 non-null  float64       
 8   Longitude              388756 non-null  float64       
dtypes: datetime64[ns](1), float64(2), object(6)
memory usage: 29.7+ MB

Loại bỏ các giá trị trùng lặp

data.duplicated().sum()

24533

data.drop_duplicates(inplace=True)
data.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 364223 entries, 95308704134 to 95316305151
Data columns (total 9 columns):
 #   Column                 Non-Null Count   Dtype         
---  ------                 --------------   -----         
 0   Incident Datetime      364223 non-null  datetime64[ns]
 1   Incident Day of Week   364223 non-null  object        
 2   Incident Category      364223 non-null  object        
 3   Resolution             364223 non-null  object        
 4   Intersection           364223 non-null  object        
 5   Police District        364223 non-null  object        
 6   Analysis Neighborhood  364223 non-null  object        
 7   Latitude               364223 non-null  float64       
 8   Longitude              364223 non-null  float64       
dtypes: datetime64[ns](1), float64(2), object(6)
memory usage: 27.8+ MB

Đặt tên lại cho các cột

data.rename(columns={
    'Incident Datetime': 'Datetime', 
    'Incident Day of Week': 'Weekday',
    'Incident Category': 'Category',
    'Police District': 'District',
    'Analysis Neighborhood': 'Neighborhood',
}, inplace=True)
data.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 364223 entries, 95308704134 to 95316305151
Data columns (total 9 columns):
 #   Column        Non-Null Count   Dtype         
---  ------        --------------   -----         
 0   Datetime      364223 non-null  datetime64[ns]
 1   Weekday       364223 non-null  object        
 2   Category      364223 non-null  object        
 3   Resolution    364223 non-null  object        
 4   Intersection  364223 non-null  object        
 5   District      364223 non-null  object        
 6   Neighborhood  364223 non-null  object        
 7   Latitude      364223 non-null  float64       
 8   Longitude     364223 non-null  float64       
dtypes: datetime64[ns](1), float64(2), object(6)
memory usage: 27.8+ MB

Tạo một biến phân loại cho thời gian

def daypart(hour):
    if   hour in ['23','00','01','02']: return 'Midnight'
    elif hour in ['03','04','05','06']: return 'Early Morning'
    elif hour in ['07','08','09','10']: return 'Morning'
    elif hour in ['11','12','13','14']: return 'Noon'
    elif hour in ['15','16','17','18']: return 'Evening'
    return 'Night'

data['Daypart'] = data['Datetime'].dt.time.apply(lambda x: daypart(str(x).split(':')[0]))
data.head()

	Datetime	Weekday	Category	Resolution	Intersection	District	Neighborhood	Latitude	Longitude	Daypart
Row ID
95308704134	2020-08-15 12:43:00	Saturday	Assault	Open or Active	GENEVA AVE \ LONDON ST	Ingleside	Excelsior	37.716039	-122.440255	Noon
95319604083	2020-08-16 03:13:00	Sunday	Assault	Open or Active	23RD ST \ ARKANSAS ST	Bayview	Potrero Hill	37.754827	-122.397729	Early Morning
95326228100	2020-08-16 03:38:00	Sunday	Malicious Mischief	Open or Active	VALENCIA ST \ 15TH ST	Mission	Mission	37.766540	-122.422044	Early Morning
95336264020	2020-08-16 13:40:00	Sunday	Non-Criminal	Open or Active	04TH ST \ MINNA ST	Southern	Financial District/South Beach	37.784044	-122.403712	Noon
95335012010	2020-08-16 16:18:00	Sunday	Weapons Offense	Cite or Arrest Adult	ORTEGA ST \ 48TH AVE	Taraval	Sunset/Parkside	37.751003	-122.507416	Evening

Phát hiện và loại bỏ outliers

plt.figure(figsize=(10, 5))

plt.subplot(1,2,1)
plt.boxplot(data['Latitude'])
plt.title('Latitude', fontsize=15)

plt.subplot(1,2,2)
plt.boxplot(data['Longitude'])
plt.title('Longitude', fontsize=15)

plt.show()

import scipy
for col in ['Latitude', 'Longitude']:
    prop = data[col]
    IQR = scipy.stats.iqr(prop)

    Q1 = np.percentile(prop, 25)
    Q3 = np.percentile(prop, 75)
    
    n_O_upper = data[prop > (Q3 + 1.5 * IQR)].shape[0]
    n_O_lower = data[prop < (Q1 - 1.5 * IQR)].shape[0]
    
    outliers_per = (n_O_upper + n_O_lower) / data.shape[0]
    print('Outliers Percentage of',  prop.name, ':', outliers_per)

Outliers Percentage of Latitude : 0.006021036562765119
Outliers Percentage of Longitude : 0.06785678004958501

=> Có thể không cần loại bỏ Outliers vì số lượng không đáng kể

Exploratory Data Analysis (EDA)

plt.figure(figsize=(10, 7))

top10_category = data['Category'].value_counts()[:10]
ax = sns.countplot(
    y = 'Category', 
    data = data, 
    order = top10_category.index
)

for rect in ax.patches:
    ax.text(
        rect.get_width() - 8000, 
        rect.get_y() + rect.get_height() / 2,
        rect.get_width(),
        color = 'white',
        weight = 'bold',
        fontsize = 11
    )
    
plt.title('Top 10 loại tội phạm ở San Francisco', fontsize=16)
plt.show()

=> Tội phạm thuộc loại Larceny Theft, cao hơn đáng kể so với bất kỳ loại tội phạm nào khác

from matplotlib.animation import FuncAnimation
from IPython.display import HTML
from datetime import datetime

df = data.copy()
df['Year'] = df['Datetime'].dt.year

year_counts = df.groupby(['Year', 'Category']).count().reset_index().pivot(
    index = 'Year', 
    columns = 'Category', 
    values = 'Daypart'
)
year_counts.fillna(0, inplace=True)

def style_axes(ax):
    ax.tick_params(labelsize=5, length=0)
    ax.grid(True, axis='x', color='white')
    ax.set_axisbelow(True)
    [spine.set_visible(False) for spine in ax.spines.values()]
    
def prepare_data(df, steps=20):
    df = df.reset_index()
    df.index = df.index * steps
    
    df_expanded = df.reindex(range(df.index[-1] + 1))
    df_expanded['Year'] = df_expanded['Year'].fillna(method='ffill')
    
    df_expanded = df_expanded.set_index('Year')
    df_rank_expanded = df_expanded.rank(axis=1, method='first')
    
    df_expanded = df_expanded.interpolate()
    df_rank_expanded = df_rank_expanded.interpolate()
    
    return df_expanded, df_rank_expanded

def init():
    ax.clear()
    style_axes(ax)

def update(i):
    ax.clear()
    for bar in ax.containers: bar.remove()
    y = df_rank_expanded.iloc[i]
    
    ax.barh(
        y = y, 
        width = df_expanded.iloc[i], 
        color = plt.cm.Dark2(range(6)), 
        tick_label = df_expanded.columns
    )
    
    for rect in ax.patches:
        ax.text(
            rect.get_width() + 500, 
            rect.get_y() + rect.get_height() / 4,
            int(rect.get_width()),
            color = 'blue',
            fontsize = 5.5
        )
    
    ax.set_ylim(min(y) - 1, max(y) + 1)
    ax.set_title(
        f'Tội phạm ở San Francisco - Năm {int(df_expanded.index[i])}', 
        fontsize=10
    )

fig = plt.Figure(figsize=(6.8, 6), dpi=144)
ax = fig.add_subplot()

df_expanded, df_rank_expanded = prepare_data(year_counts)
animation = FuncAnimation(
    fig = fig, 
    func = update, 
    init_func = init, 
    frames = len(df_expanded), 
    interval = 100, 
    repeat = False
)
HTML(animation.to_jshtml())

Once Loop Reflect

Phân tích theo Resolution

color = plt.cm.winter(np.linspace(0, 10, 20))
data['Resolution'].value_counts().plot.bar(color=color)

plt.xticks(rotation=0, fontsize=11)
plt.title('Các giải pháp cho tội phạm',fontsize=16)
plt.show()

=> Hầu hết các sự cố đều đang được tiến hành xử lý

most_commons = data[data['Category'].isin(top10_category.index)]
violent = most_commons.copy()
violent['Arrest'] = np.where(violent['Resolution'].isin(['Cite or Arrest Adult', 'Unfounded']), 0, 1)

arrest_counts = violent['Category'][violent.Arrest == 1].value_counts()[:10]
total_counts = violent['Category'].value_counts()[:10]
arrest_counts = arrest_counts / (total_counts).sort_index()
total_counts = total_counts / (total_counts).sort_index()

total_counts.plot.barh(color='crimson', label= 'Unsolved')
arrest_counts.plot.barh(color='mediumseagreen', label='Solved')

plt.legend(bbox_to_anchor=(1.05, 1))
plt.show()

plt.figure(figsize = (12, 10))

arrest = data[data['Resolution'].isin(['Cite or Arrest Adult', 'Cite or Arrest Juvenile'])]
a = arrest['Category'].value_counts().reset_index()
b = data['Category'].value_counts().reset_index()

a = a.merge(b, how='inner', on='index')
a.columns = ['Category', 'Arrests', 'Cases']

a['Arrests Percent'] = round(100 * (a['Arrests'] / a['Cases']), 2)
a.sort_values('Arrests Percent', ascending = False, inplace=True)

sns.barplot(x='Arrests Percent', y='Category', data=a)
plt.xlabel('')
plt.ylabel('')

plt.title('Điều gì sẽ khiến ta bị bắt ở San Francisco', fontsize=16)
plt.show()

=> Ma túy, vi phạm giao thông và lệnh bắt giữ là 3 lý do hàng đầu cho các vụ bắt giữ

Phân tích theo Quận

plt.figure(figsize=(10, 7))
color = plt.cm.spring(np.linspace(0, 1, 12))
ax = data['District'].value_counts().plot.bar(color=color)

for p in ax.patches:
    height = p.get_height()
    ax.text(
        x = p.get_x() + p.get_width() / 2, 
        y = height + 1000, 
        s = height, 
        ha = 'center'
    )
    
plt.xticks(rotation=0, fontsize=11)
plt.title('Quận có nhiều tội phạm nhất', fontsize=16)
plt.show()

=> Quận Central là nơi ghi nhận nhiều tội phạm nhất và quận Park là nơi ít tội phạm nhất ở San Francisco

df = pd.crosstab(data['Category'], data['District'])
color = plt.cm.Greys(np.linspace(0, 1, 10))
df.div(df.sum(1).astype(float), axis=0).plot.bar(stacked=True, color=color, figsize=(15, 7))

plt.legend(loc='upper left', bbox_to_anchor=(1, 0, 0.5, 1), fontsize=12)
plt.xlabel('')

plt.title('Quận vs Loại tội phạm', fontsize=16)
plt.show()

plt.figure(figsize = (12, 5))

a = arrest['District'].value_counts().reset_index()
b = data['District'].value_counts().reset_index()

a = a.merge(b,how = 'inner', on = 'index')
a.columns = ['District', 'Arrests','Cases']

a['Arrests Percent'] = round(100 * a['Arrests'] / a['Cases'], 2)
a.sort_values('Arrests Percent', ascending = False, inplace = True)

sns.barplot(
    x = 'District', 
    y = 'Arrests Percent',
    data = a, 
    palette = ['red'] + ['grey'] * 8 + ['blue'] + ['grey']
)

for i in range(10):
    plt.text(
        x = i, 
        y = a['Arrests Percent'] .iloc[i], 
        s = f"{a['Arrests Percent'].iloc[i]}%",
        horizontalalignment = 'center', 
        verticalalignment = 'bottom',
        weight = 'bold'
    )

plt.xticks(fontsize=11)    
plt.xlabel('')
plt.ylabel('')

plt.title('Nơi mà cảnh sát nghiêm khắc nhất - Phần trăm bị bắt', fontsize=16)
plt.show()

=> Mặc dù có ít tội phạm hơn được ghi nhận ở quận Tenderlion nhưng tỷ lệ bắt giữ ở đây khá cao 36,09%

Phân tích theo Thời gian

del df
df = data.copy()
df = df.groupby(df['Datetime'].dt.date).count().iloc[:, 0]

color_palette = sns.color_palette()
sns.kdeplot(data=df, shade=True)

plt.axvline(
    x = df.median(), 
    ymax = 0.95, 
    linestyle = '--', 
    color = color_palette[1]
)

plt.annotate(
    f'Median: {df.median()}',
    xy = (df.median(), 0.004),
    xytext = (200, 0.005),
    arrowprops = dict(
        arrowstyle = '->', 
        color = color_palette[1], 
        shrinkB = 10
    )
)

plt.xlabel('Incidents')
plt.ylabel('Density')

plt.title('Phân phối số lượng sự cố mỗi ngày', fontdict={'fontsize': 16})
plt.show()

cat_per_week_common = pd.crosstab(most_commons['Category'], most_commons['Weekday']) 
cat_per_week_common = cat_per_week_common.div(cat_per_week_common.sum(axis=1), axis=0)
cat_per_week_common = cat_per_week_common[['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday','Sunday']]

ax = sns.heatmap(cat_per_week_common, cmap="BuPu", linewidths=0.5)
plt.xticks(fontsize=11, rotation=45, ha='right')
plt.yticks(fontsize=11)

plt.xlabel('')
plt.ylabel('')
plt.show()

data['Weekday'].value_counts().plot.pie(
    explode = [0.1, 0, 0, 0, 0, 0, 0], 
    textprops = {'color': 'w', 'fontsize': 15, 'weight': 'bold'},
    autopct = "%.1f%%"
)

plt.legend(loc='center left', bbox_to_anchor=(1, 0, 0.5, 1), fontsize=12)
plt.axis('off')

plt.title('Số lượng tội phạm theo thứ', fontsize=16)
plt.show()

=> Thứ 6 là ngày mà tội phạm được ghi nhận nhiều nhất, tiếp theo là Thứ 4. Chủ nhật là ngày ít tội phạm nhất

del df
df = data.copy()
df['Date'] = df['Datetime'].dt.date
df['Hour'] = df['Datetime'].dt.hour

df = df.groupby(['Hour', 'Date', 'Category'], as_index=False).count().iloc[:, :4]
df.rename(columns={'Datetime': 'Incidents'}, inplace=True)

df = df.groupby(['Hour', 'Category'], as_index=False).mean()
df = df[df['Category'].isin(top10_category[4:].index)]

fig, ax = plt.subplots(figsize=(12, 5))
ax = sns.lineplot(x='Hour', y='Incidents', data=df, hue='Category')
ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15), ncol=6)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

plt.suptitle('Số sự cố trung bình mỗi giờ')
plt.show()

plt.figure(figsize=(15, 7))
color = plt.cm.twilight(np.linspace(0, 4, 100))
data['Datetime'].dt.time.value_counts().head(24).sort_index().plot.bar(color=color)

plt.xticks(rotation=45, fontsize=11)
plt.title('Phân bố tội phạm trong ngày', fontsize=16)
plt.show()

df = data['Daypart'].value_counts().reset_index()
df_normalize = data['Daypart'].value_counts(normalize=True)
sns.barplot(x=df['index'], y='Daypart', data=df)

for i in range(6):
     plt.text(
         x = i, 
         y = 5000, 
         s = f'{round(100 * df_normalize[i], 2)}%', 
         horizontalalignment = 'center', 
         color = 'white',
         weight = 'bold'
     )
        
plt.title('Số lượng tội phạm theo các buổi trong ngày', fontsize=16)
plt.show()

=> Hãy cẩn thận vào buổi tối. Gần 1/4 số vụ tội phạm xảy ra vào buổi tối

plt.figure(figsize = (15, 7))

pivot_table = pd.pivot_table(
    columns = data['Daypart'] ,
    index = 'Weekday',
    values = 'Daypart' , 
    aggfunc = 'count', 
    data = data
)

pivot_table = pivot_table.reindex(
    index = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday'], 
    columns = ['Early Morning', 'Morning', 'Noon', 'Evening', 'Night', 'Midnight']
)

sns.heatmap(pivot_table, cmap='Reds', annot=True, fmt='d')
plt.xlabel('')
plt.ylabel('')
plt.show()

=> Buổi tối thường là lúc có nhiều vụ phạm tội, tối Thứ 6 đặc biệt tồi tệ và tối Chủ nhật tương đối tốt hơn

sns.countplot(data['Datetime'].dt.month, palette='autumn')
plt.xticks(fontsize=11)
plt.xlabel('')

plt.title('Tội phạm theo từng tháng', fontsize=16)
plt.show()

=> Tháng 1 là tháng có nhiều tội phạm nhất

plt.figure(figsize = (15, 7))

df = data[['Datetime']]
df['Day'] = df['Datetime'].dt.day
df['Month'] = df['Datetime'].dt.month

pivot_table = pd.pivot_table(
    columns = df.Day, 
    values = 'Day', 
    index = 'Month', 
    aggfunc = 'count', 
    data = df
)
sns.heatmap(pivot_table, cmap='Reds')

plt.xticks(fontsize=12)
plt.yticks(
    np.arange(0.5, 12.5),
    ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sept', 'Oct', 'Nov', 'Dec'],
    fontsize = 12
)

plt.xlabel('')
plt.ylabel('')
plt.show()

=> Ngày 1/1 dường như là ngày có số tội phạm được ghi nhận cao nhất và ngày 25/12 là ngày có số tội phạm được ghi nhận ít nhất

Phân tích Time Series

first_date = data['Datetime'].min()
last_date = data['Datetime'].max()

print('Ngày đầu tiên:', first_date)
print('Ngày cuối cùng:', last_date)
print(f'Số ngày tổng cộng: {(last_date - first_date).days + 1}')

Ngày đầu tiên: 2018-01-01 00:00:00
Ngày cuối cùng: 2020-12-04 22:19:00
Số ngày tổng cộng: 1069

time_span = pd.date_range(first_date, last_date)
def to_time_series(df):
    date = df['Datetime'].value_counts().index
    ts = pd.DataFrame(
        data = df['Datetime'].value_counts().values[np.argsort(date)], 
        index = sorted(date), 
        columns = ['count']
    )
    return ts.reindex(time_span, fill_value=0) 

all_counts = to_time_series(data)
all_counts[all_counts['count'] == 0]

	count

=> Không có ngày nào là bình yên (không có bất kỳ sự cố nào xảy ra)

plt.figure(figsize=(16, 7))
all_counts['count'].plot()

plt.xlabel('Thới gian')
plt.ylabel('Số lượng tội phạm')
plt.show()

=> Liệu các thay đổi này có đúng cho các loại tội phạm riêng lẻ ? Ta sẽ thử với Larceny Theft

df_theft = data[data['Category'] == 'Larceny Theft']
theft_counts = to_time_series(df_theft)
theft_counts[theft_counts['count'] == 0]

	count
2018-01-03	0
2018-01-11	0
2018-02-05	0
2018-02-08	0
2018-02-22	0
...	...
2020-10-24	0
2020-10-28	0
2020-11-10	0
2020-11-28	0
2020-11-29	0

65 rows × 1 columns

plt.figure(figsize=(16, 7))
theft_counts['count'].plot()

plt.xlabel('Thới gian')
plt.ylabel('Số lượng tội phạm Larceny Theft')
plt.show()

=> Nó thực sự khác nhau nhiều

y = theft_counts.resample('MS').sum()
y.tail()

	count
2020-08-01	78
2020-09-01	77
2020-10-01	72
2020-11-01	57
2020-12-01	5

Do dữ liệu chỉ được ghi cho đến ngày 4 tháng 12 năm 2020, nên ta sẽ không dùng dữ liệu của tháng 12 năm 2020 trở đi

fig = plt.figure(figsize=(20,6))
y = y[:-1]
y['count'].plot()

plt.xlabel('Thời gian')
plt.ylabel('Số lượng Larceny Theft mỗi tháng')
plt.show()

=> Ta có thể thấy rằng có sự giảm mạnh từ tháng 7 có thể là do ảnh hưởng của Covid-19

Ở đây, những danh sách tội phạm này có thể được khám phá và mô hình hóa thêm nếu có thêm dữ kiện của nhiều năm trước:

• Ta có thể sử dụng mô hình ARIMA để huấn luyện với dữ liệu Larceny Theft.
• Để giảm tính ngẫu nhiên của dữ liệu, ta có thể phân chia chuỗi thời gian thành các khoảng thời gian hàng tháng. Chiến lược này cũng tốt nếu có sự thay đổi theo thời vụ trong dữ liệu, 12 tháng bằng 1 năm. Vì 1 năm không phải là số nguyên của tuần nên việc sử dụng khoảng thời gian hàng tuần cho mô hình ARIMA theo thời vụ không có xu hướng cho kết quả tốt.
• Loại bỏ thành phần thời vụ khỏi chuỗi thời gian có thể dẫn đến mối quan hệ rõ ràng hơn giữa các biến đầu vào và đầu ra

Phân tích theo Mô tả

data_with_description = raw_data[raw_data.index.isin(data.index)]
data_with_description.rename(columns={
    'Incident Datetime': 'Datetime', 
    'Incident Day of Week': 'Weekday',
    'Incident Category': 'Category',
    'Incident Description': 'Description',
    'Police District': 'District',
    'Analysis Neighborhood': 'Neighborhood',
}, inplace=True)

from wordcloud import WordCloud
description = data_with_description['Description']
wc = WordCloud(background_color='#e9eaf1', width=1000, height=500).generate(str(description))

plt.imshow(wc)
plt.axis('off')

plt.title('Mô tả cho các tội', fontsize=20)
plt.show()

import plotly.express as px
discat = data_with_description.groupby(['District', 'Category'])['Description'].count().reset_index()
fig = px.sunburst(
    discat.rename(columns={'Description': 'Count'}), 
    path = ['District', 'Category'], 
    values = 'Count', 
    color = 'Count'
)
fig.show()

Phân tích theo Phân bố không gian (Spatial Distribution)

plt.figure(figsize=(10, 7))
ax = sns.countplot(
    y = 'Intersection', 
    data = data, 
    order = data['Intersection'].value_counts()[:15].index
)

for rect in ax.patches:
    ax.text(
        rect.get_width() - 150, 
        rect.get_y() + rect.get_height() / 2,
        rect.get_width(),
        color = 'white',
        weight = 'bold', 
        fontsize = 11
    )
    
plt.title('Top 15 vùng phạm tội', fontsize=16)
plt.show()

del df
df = data.copy()

mean_lat = np.mean(df.Latitude)
mean_lon = np.mean(df.Longitude)

sw = df[['Latitude', 'Longitude']].min().values.tolist()
ne = df[['Latitude', 'Longitude']].max().values.tolist()

df['Coordinates'] = df['Longitude'].astype(str) + ', ' + df['Latitude'].astype(str)
df.Coordinates.value_counts()[:5]

-122.40733700000001, 37.78456014    2490
-122.4036355, 37.77516081           2298
-122.4080362, 37.78640961           2160
-122.419669, 37.76505134            1681
-122.41259529999999, 37.78393258    1504
Name: Coordinates, dtype: int64

=> Cặp tọa độ phổ biến nhất (-122.40733700000001, 37.78456014) là vị trí của Market St

plt.title('Phân bố tội phạm')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.scatter(data.Longitude, data.Latitude, s=0.5, c='r')
plt.show()

=> Ta có thể thấy có nhiều tội phạm hơn ở phần đông bắc của San Francisco. Câu hỏi đặt ra rằng liệu điều này có luôn như vậy trong suốt các năm qua không ?

Ta sẽ xem liệu phân phối có bất kỳ thay đổi đáng kể nào không:

• Để giảm thời gian tính toán, ở đây sẽ chỉ vẽ dữ liệu của tháng 1 (tháng có nhiều tội phạm nhất) hàng năm.
• Để hiểu rõ hơn về sự phân bố không gian, lần này biểu đồ sẽ được phủ các đường bao lên các ô phân tán.

plt.figure(figsize = (16, 5))
for index, year in enumerate(range(2018, 2021)): 
    df_jan = data[(data['Datetime'] < f'{year}-02') & (data['Datetime'] > f'{year - 1}-12-31')]
    
    plt.subplot(1, 3, index + 1)
    plt.scatter(df_jan.Longitude, df_jan.Latitude, s=0.5, c='r')
    sns.kdeplot(df_jan.Longitude, df_jan.Latitude)
    
    plt.xticks(rotation=45)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    
    plt.title(f'Phân bố tội phạm vào tháng 1 năm {year}')
plt.show()

=> Có vẻ như phần đông bắc của San Francisco luôn là khu vực nguy hiểm nhất

Trực quan hóa không gian địa lý

import geopandas as gpd
from shapely.geometry import Point

gdf = data.copy()
gdf['Coordinates'] = list(zip(gdf.Longitude, gdf.Latitude))
gdf.Coordinates = gdf.Coordinates.apply(Point)

gdf = gpd.GeoDataFrame(gdf, geometry='Coordinates', crs={'init': 'epsg:4326'})
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

ax = world.plot(color='white', edgecolor='black')
gdf.plot(ax=ax, color='red')
plt.show()

district_labels = {
    'Southern': '37.774432, -122.401121',
    'Bayview': '37.734332, -122.389920',
    'Mission': '37.756478, -122.423663',
    'Northern': '37.787740, -122.430300',
    'Tenderloin': '37.781980, -122.412981',
    'Central': '37.796200, -122.409293',
    'Park': '37.765352, -122.449282',
    'Richmond': '37.776204, -122.483285',
    'Ingleside': '37.726817, -122.437207',
    'Taraval': '37.737775, -122.484375'
}

table = data['District'].value_counts().reindex(district_labels.keys())
table = table.reset_index().rename(
    {'index': 'District', 'District': 'Count'}, 
    axis='columns'
)

table['District'] = table['District'].str.upper()
table

	District	Count
0	SOUTHERN	44896
1	BAYVIEW	32120
2	MISSION	48846
3	NORTHERN	50818
4	TENDERLOIN	35127
5	CENTRAL	55065
6	PARK	17638
7	RICHMOND	21415
8	INGLESIDE	28329
9	TARAVAL	25815

import urllib.request
import shutil
import zipfile

url = 'https://data.sfgov.org/api/geospatial/wkhw-cjsf?method=export&format=Shapefile'
with urllib.request.urlopen(url) as response, open('pd_data.zip', 'wb') as out_file:
    shutil.copyfileobj(response, out_file)
    
with zipfile.ZipFile('pd_data.zip', 'r') as zip_ref:
    zip_ref.extractall('pd_data')

import os
import re
    
for filename in os.listdir('./pd_data/'):
    if re.match(".+\.shp", filename):
        districts = gpd.read_file(f'./pd_data/{filename}')
        break

districts.crs={'init': 'epsg:3857'}
districts = districts.merge(
    table.set_index(['District']),
    how = 'inner',
    left_on = 'district',
    right_index = True,
    suffixes = ('_x', '_y')
)
districts

	company	district	shape_area	shape_le_1	shape_leng	geometry	Count
0	B	SOUTHERN	9.134414e+07	100231.353916	87550.275142	MULTIPOLYGON (((-122.39186 37.79425, -122.3917...	44896
1	C	BAYVIEW	2.013846e+08	144143.480351	163013.798332	POLYGON ((-122.38098 37.76480, -122.38103 37.7...	32120
2	D	MISSION	8.062384e+07	40518.834235	40152.783389	POLYGON ((-122.40954 37.76932, -122.40862 37.7...	48846
3	E	NORTHERN	8.278169e+07	50608.310321	56493.858208	POLYGON ((-122.43379 37.80793, -122.43375 37.8...	50818
4	J	TENDERLOIN	1.107215e+07	18796.784185	12424.268969	POLYGON ((-122.40217 37.78626, -122.41718 37.7...	35127
5	A	CENTRAL	5.595027e+07	67686.522865	64025.129073	POLYGON ((-122.42612 37.80684, -122.42612 37.8...	55065
6	F	PARK	8.487896e+07	50328.913294	46307.776968	POLYGON ((-122.43956 37.78314, -122.43832 37.7...	17638
7	G	RICHMOND	1.379640e+08	75188.628361	69991.465355	POLYGON ((-122.44127 37.79149, -122.44060 37.7...	21415
8	H	INGLESIDE	1.935805e+08	74474.181164	74737.936295	POLYGON ((-122.40450 37.74858, -122.40407 37.7...	28329
9	I	TARAVAL	2.846767e+08	73470.424000	75350.217521	POLYGON ((-122.49842 37.70810, -122.49842 37.7...	25815

import contextily as ctx
districts['incidents_per_day'] = districts.Count / data.groupby('Datetime').count().shape[0]

fig, ax = plt.subplots(figsize=(12, 7))
districts.plot(
    column = 'incidents_per_day',
    cmap = 'Reds',
    alpha = 0.6,
    edgecolor = 'r',
    linestyle = '-',
    linewidth = 1,
    legend = True,
    ax = ax)

for index in districts.index:
    name = districts.loc[index].district
    geometry = districts.loc[index].geometry
    
    plt.annotate(
        name,
        (geometry.centroid.x, geometry.centroid.y),
        color = '#353535',
        fontsize = 'large',
        fontweight = 'heavy',
        horizontalalignment = 'center'
    )
plt.show()

import geoplot as gplt
crimes = data['Category'].unique().tolist()

sf_land = districts.unary_union
sf_land = gpd.GeoDataFrame(gpd.GeoSeries(sf_land), crs={'init':'epsg:4326'})
sf_land = sf_land.rename(columns={0: 'geometry'}).set_geometry('geometry')

fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(15, 7))
for i, crime in enumerate(np.random.choice(crimes, size=6, replace=False)):
    ax = fig.add_subplot(2, 3, i + 1)
    gplt.kdeplot(
        gdf.loc[gdf['Category'] == crime],
        shade = True,
        shade_lowest = False,
        clip = sf_land.geometry,
        cmap = 'Reds',
        ax = ax)
    gplt.polyplot(sf_land, ax=ax)
    ax.set_title(crime) 

fig.tight_layout(rect=[0, 0, 0, 0])
ax.set_axis_off()

plt.title('Mật độ địa lý của các tội phạm khác nhau')
plt.show()

import folium
gjson = r'https://cocl.us/sanfran_geojson'

figure = folium.Figure(width=900, height=600)
sf_map = folium.Map(location=[37.77, -122.42], zoom_start=12)
sf_map.choropleth(
    geo_data = gjson,
    data = table,
    columns = ['District', 'Count'],
    key_on = 'feature.properties.DISTRICT',
    fill_color = 'YlOrRd', 
    fill_opacity = 0.7, 
    line_opacity = 0.2,
    legend_name = 'Tỉ lệ tội phạm ở San Francisco'
)
    
figure.add_child(sf_map)
figure

%%html
<iframe 
    src="https://data.sfgov.org/dataset/Map-of-Police-Department-Incident-Reports-2018-to-/jq29-s5wp/embed?width=950&height=600" 
    width="950" 
    height="600" 
    style="border:0; padding: 0; margin: 0;"
></iframe>

Step 3 : Model selection & Trainning

Step 4: Evaluate & test the model results & performance

Step 5: Deploy the model